;; LTDL_LIBRARY_PATH=/opt/local/lib
;; GUILE_LOAD_PATH=.

;;(import (srfi srfi-1))

(import (rnrs) 
        (grid-gen common)
)

(module-on fortran77)

;; (module-on real-8)
;; (module-on OpenCL)
;; (module-on OpenCL-CPU)

(module-on HD)

;;(module-on Lax-Friedrichs)
(module-on Roe)
(module-on Einfeldt)
;;(module-on Harten2)
(module-on Osher)

(module-on grid)

(module-on OpenMP #t)
;; (module-on MPI 4)
(module-on maxima)
(module-on OpenDX)

(define tim (comp-float "tim" 0.0))
(define tau (comp-float "tau"))

(define grid-dims '(100 100 200))
;;(define grid-dims '(1000 2000))

(define pi 3.1415926)

(def-float xmin -1.0)
(def-float xmax 1.0)

(def-float ymin -1.0)
(def-float ymax 1.0)

(def-float zmin 0.0)
(def-float zmax 4.0)

(def-float tmax 0.8)

(def-float sr 0.1)

(def-float rho0 1.0)
(def-float p0     1.0)
(def-float u0     0.0)
(def-float v0     0.0)
(def-float w0     0.0)

(def-float rho1 10.0)
(def-float p1     1.0)
(def-float u1     0.0)
(def-float v1     0.0)
(def-float w1     5.0)

(def-float cour  0.1)

(define gr (grid:cart grid-dims))

(define gr-cond `(< (+ (** ,(gr 'x) 2.0) (** ,(gr 'y) 2.0)) (** ,sr 2.0)))

(gr 'add-bnd 'z-
   (beg
     (set-val (gr 'r)   `(if ,gr-cond ,rho1 ,rho0))
     (set-val (gr 'p)   `(if ,gr-cond ,p1   ,p0))
     (set-val (gr 'u)   `(if ,gr-cond ,u1   ,u0))
     (set-val (gr 'v)   `(if ,gr-cond ,v1   ,v0))
     (set-val (gr 'w)   `(if ,gr-cond ,w1   ,w0))
    )
)

(gr 'add-corr 'smooth)
(gr 'add-corr 'err)

(define grdf0 (lambda (x y z)
  (list `(+ ,xmin (* (- ,xmax ,xmin) ,x))
        `(+ ,ymin (* (- ,ymax ,ymin) ,y))
        `(+ ,zmin (* (- ,zmax ,zmin) ,z))
  )
))

(def-float G 3.0)

(define mypot (lambda (x y z)
  `(* ,x ,G)
))

;;(set! PHYS:Potential mypot)

(gr 'set-geom grdf0)

(code-gen "stream3d"
  (beg
    (gr 'init-grid)
    (gr 'init rho0 p0 u0 v0 w0)

    (do-while (and (< tim tmax) (= (gr 'err) 0))

      (gr 'update-bnd)
      
      (if (= (gr 'err) 0)
        (beg
          (gr 'CFL cour)
          (if (> (+ (gr 'tau) tim) tmax)
              (set-val (gr 'tau) (- tmax tim))
          )
        )
      )
      (main-p
        (do-every 1 'sec
          (print "t=" tim "tau=" (gr 'tau))
        )
      )
      (gr 'calc)
      (if (= (gr 'err) 0)
          (set-val tim (+ tim (gr 'tau)))
          (beg
            (set-val (gr 'tau) (* 0.5 (gr 'tau)))
            (print "new Tau=" (gr 'tau))
          )
      )
    )
    (OpenDX:save gr "stream3d.dx")
  )
)

(display "\n")
